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Abstract 

We present here a general iterative formula which gives a (formal) 
series expansion for the time autocorrelation of smooth dynamical vari- 
ables, for all Hamiltonian systems endowed with an invariant measure. 
We add some criteria, theoretical in nature, which enable one to de- 
_C . cide whether the decay of the correlations is exponentially fast or not. 

One of these criteria is implemented numerically for the case of the 
Fermi-Pasta-Ulam system, and we find indications which might sug- 
gest a sub-exponential decay of the time autocorrelation of a relevant 
dynamical variable. 

1 Introduction 

It is well known that one of the most important indicators of the chaotic 
behaviour of a dynamical system is the rate of decay to zero of the time 
autocorrelation of dynamical variables. A vast literature exists on this sub- 
ject, mainly addressed to systems of some special class (for example, Anosov 
systems, see [T], or systems with few degrees of freedom, such as billiards, 
see [2]). In this paper we study the decay of correlations in the general frame 
of the dynamical theory of Hamiltonian systems, with an invariant measure. 
Such an approach has already provided interesting results, as it allows to 
obtain results valid in the thermodynamic limit (see |3j SI E] ) - 

Here we provide a power series expansion (with respect to time) of the 
time autocorrelation of a smooth function / (Theorem [2l Section [2]). It turns 
out that the coefficients of such a series are essentially the variances of Vjjf, 
where Lh is the Poisson bracket operator relative to the Hamiltonian H. We 
then establish a necessary and sufficient condition for its time-decay not to be 
exponential (Corollary Q] Section [3]). We finally give a numerical application 
to the case of the Fermi-Pasta-Ulam system (FPU). By truncating the series 
up to order 12, for systems up to 364 degrees of freedom, we approximate 
the poles of the Laplace transform of the time autocorrelation of a suitably 
chosen dynamical variable. Some comments are made on the possibility 
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that such approximations give indications for a sub-exponential decay of the 
autocorrelation. 

In Section [2] the relevant notions on time correlation functions are re- 
called. In the same section we give a theorem which establishes a kind of 
continuity (with respect to a suitable distance) of time correlations in the 
L 2 space of the dynamical variables, and also Theorem [2] for the series ex- 
pansion. 

Section [3] is devoted to a study of some general properties of the series 
expansion. Here an interesting link with the Stieltjes moment problem is 
pointed out (see Theorem [3]), and a necessary and sufficient criterion for the 
time-decay of the autocorrelation to be exponential is given as Corollary [TJ 
A closer examination of this problem is given in Appendix lAl 

Then, the results of some numerical computations on the time autocor- 
relations of a significant variable in a FPU chain are reported in Section 
Finally, in Section [5] some comments on the usefulness of the present method 
are also given. Two more Appendices complete the paper. 

2 Time correlations in their "natural" space 

We recall here some standard concepts within the measure theoretic approach 
to dynamical systems. We consider a Hamiltonian system on a phase space 
A4 endowed with a probability measure fx invariant with respect to the time 
flow 4> t induced by the Hamiltonian H. We also consider the time evolution 
operator Ut acting on the space of the square integrable functions from Ai 
to M, i.e., L 2 (n,A4). The operator Ut maps / to Utf = f o It is known 
after Koopman (see [6]) that Ut defines a one-parameter group of unitary 
operators, namely operators preserving the norm || • || in L 2 (/j,,M). 

We are interested in the time autocorrelation of a dynamical variable 
/ G L 2 (fi,Ai), which is defined as 

C f (t) d ^ f (ft f) - (f) 2 = ((ft - (/)) (/ - (/))) , (1) 

where ft is a shorthand for Utf, and (•) denotes mean value with respect to 
the probability measure \i. 

It will also be useful to consider the time correlation between two dy- 
namical variables / and g, defined as 

C f , 9 (t) d ^ f (f t g) - (f)(g) = ((ft - (/» (g - (g))) . 

It is well known from probability theory that the concepts of variance 
and correlation acquire a geometrical meaning if the covariance between two 
random variables / and g is used as a scalar product. The covariance is 
defined as 

cov(/, g) d ^ f ((/-(/}) (g-(g))) = C fig (0) , 
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This choice leads to take as nor of a dynamical variable / its standard 
deviation 07, defined by 

^ d 4 f cov(/,/) = ||/|| 2 -(/) 2 = ||/-(/)|| 2 , 

cr 2 being the variance of /. 

The notion that some relation exists between two random variables / 
and g is made quantitative through the correlation coefficient r(/, g) defined 
as 

r(/,3) = cov{f,g)/a f a g . 

Using the Schwarz inequality it is easily seen in our case that one has 
|cov(/, g)\ < &fCTg, so that |r(/, g)\ < 1. Two variables are orthogonal 
in this metric if they are uncorrelated, i.e., if r(/, g) = 0, and collinear if 
r(/,<?)=±l. 

We also emphasise that, in view of the above definitions, the time auto- 
correlation of a function / is nothing but the covariance between Utf and /. 
Therefore one immediately has 

0/(0) = a) and |C/(t)| < a) . (2) 

Thus, the variance of a dynamical variable is the natural scale of its time 
autocorrelation. 

We come now to state two useful theorems. The first one points out an 
interesting property of the time correlation. The mere fact that two dynam- 
ical variables are strongly correlated (i.e., sufficiently collinear in the scalar- 
product given by the covariance) entails a similar behaviour of their time 
autocorrelations. Equivalently, we can say that in the neighbourhood of 
each dynamical variable /, i.e., in the set of all dynamical variables strongly 
correlated with /, the behaviour of the time autocorrelations is determined 
by that of C f (t). 

Theorem 1 Let \x be a probability measure on the phase space Ai, invariant 
for the flow generated by H , and let f and g be dynamical variables belonging 
to L 2 (fj,,Ai) such that one has |r(/, g)\ > 1 — e 2 /2 for some e > 0. Then 

there exists a multiple f = af of f , with a = f (sign(x(f ,g))a g / a f), such that 

\C f ~ g (t)-C f (t)\<ea 2 g (3) 

and 

|C^)-C / (t)|<(e 2 + 2e)a 2 . (4) 

1 We can adopt this quantity as a norm in strict sense only if all dynamical variables 
which differ by a constant are identified. For a function with zero mean the covariance 
cov(/, /) actually coincides with the usual L 2 norm. In this paper we shall use both norms, 
keeping the usual symbol || • || for the L 2 norm. 
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Proof. Both inequalities come from the remark that 

a g _- f = (T 2 g + crj - 2cov(/,s) = 2a 2 - 2a 2 x(f,g) < e 2 a 2 , 

which is due to the identities = a g and r(f,g) = |r(/, ^)|. In fact, ([3]) 
hence follows by noting that 

\ C f,g^> ~ C /(*)| = \ cow Ut,g) ~ cov(/ t ,/)| = |cov(/ t ,5- /)| < <r g cr g _f , 
while, in a similar way, @ comes from the following relations 

\C g (t) - C f (t)\ = \C g _ f (t) + C gJ (t) + C lg (t) - 2C / (t)| 

K < J • H% /• 

Here, in the second line, use is made of the previous inequality and of the 
identity C g ^(i) = C s (— t), due to the invariance of the measure. 

Q.E.D. 

A central role in this paper will be played by the formal power series 
expansion of Cf(t) with respect to time, which is given in the following 
theorem. Use will be made of the definition of the fc-th order Lie derivative 
of /, for k > 1, namely d = [/( fc-1 ), il], where = f / and [•, •] denotes 
the Poisson brackets. 

Theorem 2 Let \i be a probability measure on the phase space A4, invariant 
for the flow generated by H , and let n > 0. Then, for any dynamical variable 
f G L 2 (n,M) such that f^ £ L 2 (^,M) for all k <n, one has 



(k) 



,2 t 2k 



+(-ir +i fdh... r n_i di 2n |i/t)-/wi 

Jo Jo 



(5) 



Proof. We use the simple chain of identities 

\\f t -f\\ 2 = 2a 2 -2C f (t) = 2\\f\\ 2 -2(f t f) , (6) 
which hold true for any invariant measure, and the remark that 

ft~f= fut-sfds, (7) 
Jo 

where / = [/, H] . One can check equations ([6]) by writing down the square at 
the l.h.s., while (J7|) comes H from the variation of constants formula applied 
to f t - f. 



2 See |2| for a more detailed proof. 
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We go on by writing the l.h.s. of ([6]) as 

ll/t-/|| 2 = fdfi ! dt' ! df(U t „ t ,f)(U t „ t „f) 



o jo 

t ft 



dfi [ dt' [ dt" (8) 
J o J o 

= * f& f (\\f§ - ft\ ( 9 ) 

where, in the second line, the invariance of the measure with respect to the 
time evolution is used, so that the equality proceeds from a simple change of 
coordinates, and in the last line the same argument is used, as well as relation 
([6|) applied to the variable /. The order of the integrals can be exchanged 
if ll/H is finite, because the integral Q is absolutely convergent in this case, 
in virtue of Schwarz inequality. In @ we repeat the same calculation for 
Wft" — /|| 2 j using / = [f,H]. This gives us the relation 

2t 2 „,.™„2t 4 



II/* -/II 2 = 2(||/|| 2 |-||/( 2 ) 



4! . 

pi pt\ rt'2 rts 

+ dh dt 2 / dt 3 / dt 4 ||A (2) -/ (2) || 2 , 

Jo Jo Jo Jo 

where /( 2 ) = [f,H]. The same procedure can be iterated up to k = n by 

recursively defining ^ [/( fe_1 ) , ff] , starting with /(°) = f /, and using 
the hypothesis that £ L 2 (p,A4). This completes the proof. 

Q.E.D. 

Remark 1. Theorem [2] suggests that, at least in a formal way, one can 
express the time autocorrelation of a dynamical variable / as a power series 
with respect to time 

+°° + 2k 

C/(t)=4 + EH fe ||/ (fe) ll (k)T- ( 10 ) 
To simplify the notations, we define 

co^aj, c n d ^ f ||/W|| 2 forn>0. (11) 

We point out that all these coefficients are positive. We could have got 
to formula (fTU|) also by expressing ft as a time power series through its 
Lie derivatives and by integrating (ft f) by parts, taking into account the 
observation that the mean value of any function which can be written as 
[g,H] vanishes^, since {[g,H]) = ^(g). 



3 For the same reason, <?.(*:) can replace in all the previous formulae. We will 

keep the given notation, because this way it is more straightforward to understand how 
to compute the coefficients of the series. 
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Remark 2. We emphasise that the term in the second line in equation 
(|5|), i.e., the remainder, turns out to be positive or negative according to 
n being odd or even, respectively. Therefore, the truncations of series (|10p 
provide upper bounds for the time autocorrelation of / if truncated at an 
odd n, and lower bounds if truncated at an even n. Such bounds provide 
some information on the behaviour of the time autocorrelation of / for finite 
times. The simplest example is that of the first order truncation, which 
proved to be very helpful for the study of relaxation times in a Hamiltonian 
setting (see [31110). 

As a final remark, we point out that some a priori restrictions on the 
norms of the successive derivatives come from the Schwarz inequality Q. 
Indeed, one has the following Proposition, whose proof is deferred to Ap- 
pendix [Bj 

Proposition 1 The coefficients c n defined by Ml\) are such that the polyno- 
mials Pn{y), Q l n {y) defined as 

2n 2n 

f*M d ^ 2a l + , QUy) = £(-M + i/ , (12) 

k=l k=0 

with 

a l n ^ c n+l /(2n)\ , (13) 
are positive semi-definite, for any n,l > 0. 

The statement of Proposition [1] entails some restrictions on the coefficients 
c n , since P^(y) and Q l n (y) are two-indexed sequences of positive semi-definite 
polynomials. Unfortunately, such relations can be expressed only in a quite 
complicated way, through the Jacobi-Borchardt theorem (see, for example, 
[7J), which gives some relations among the roots of the polynomials, and 
Newton's identities on the relations between the coefficients of a polynomial 
and its roots. 

3 Asymptotic behaviour of the correlations 

In many cases one is interested in determining whether the time correlation 
tends to zero as t — > oo and, if this is the case, which is the law controlling 
the decay. Indeed, it is commonly stated that if the correlations decay ex- 
ponentially fast (for a suitable wide class of functions) then the system is 
chaotic. The decay is linked to the analyticity property of the Fourier trans- 
form of the correlations (see for example |8j): in particular it is well known 
that the Fourier transform is analytic in a strip of width 1/r if and only if the 
correlation decays at least as exp(— t/r). We recall that the correlations can 



always be expressed as the Fourier transform of an even positive measure, 
the so called "spectral measure", namely as 

i r+oo 

Cf(t) = — J e**da( U ) , (14) 

where a is a positive Borel measure on M + (see, for example, |llj). 

On the other hand, it is not so easy to get information about the spec- 
trum: we devise to get it using the Laplace transform F(s) of the time 
autocorrelation, namely, 

F(s) = [ +0 ° e~ st C f (t) for s £ C . (15) 
Jo 

It is well known that the Laplace transform of a function which decays faster 
than exp(— t/r) is analytic in the half-plane Res > — 1/r, so that, again, the 
control on the decay of correlations in linked to the analyticity properties of 
F(s), as in the case of the Fourier transform. The great advantage of the 
latter approach lays in the fact that the Laplace transform of the correlation 
is essentially the Stieltjes transform of the spectral measure. This is stated 
as Theorem [3] below. 

The point is that, as first shown by Stieltjes (see |9j), if one knows the 
asymptotic expansion in powers of 1/s of a functions F(s) which is the 
Stieltjes transform of a positive measure, then one can recover F(s) itself 
through a convergent scheme of continued fractions expansions, defined in 
terms of the coefficients of the asymptotic expansion. In other words, one 
can construct a rational approximation 

) ~ Qn(s) ' 

which converges uniformly in Re s > 0, as n — > +oo, to F(s). For the case of 
the time autocorrelation Cf(t), the asymptotic expansion is obtained just by 
integrating the expansion (|10p term by term. In fact, the Laplace transform 
(|15p can be formally written as 

r +oo /+°o \ +oo 

m-J, (B-r^J^^B-)^. (-) 

This expansion would give a convergent expansion, provided one can ex- 
change the sum of the series with the integration over time. This is in 
general not permissible, in view of the fact that the expansion in power se- 
ries of time converges (in general) only in a circle of finite radius and not up 
to infinity. However, it is easy to show that series (|16p nevertheless provides 
the asymptotic expansion of F(s) we are looking for. One has to use formula 
([5]) and check that the remainder grows at most as a power of time. 
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We are now able to state the main result of this section, namely, the 
following theorem, together with a relevant corollary on the decay of the 
correlations. 

Theorem 3 Let C/(i) be analytic in t about the origin and continuous for 
any t E R. Then the following statements hold: 

1. in the half-plane Res > the Laplace transform F(s) of C fit) is an- 
alytic and one has 



n*) = s -f 

7T Jo 



da(uj) 



. , > • (17) 
s z + u z 



where a{uj) is the positive Borel measure such that (|14p holds; 

2. the continued fraction which approximates the asymptotic series |7? 
converges to -F(s) for Res > and the Borel positive measure §{u) 

defined on IR + by $(it) = f {l/^T)a(^/u) solves the Stieltjes moment 
problem for the coefficients c n defined by Ml]) - Moreover, a(u)/TT solves 
the symmetric Hamburger moment problem with moments c' 2n = c 



C 2n+1 — 0- 



Corollary 1 Cf(t) decays exponentially fast as t goes to infinity if and only 
if the coefficients c n are such that a'{oj) exists and y/u^'{u) is analytic. 

Remark 3. We will state our results in terms of the regularity properties of 
a or $ according to convenience. We prefer the study of the properties of 
<£, when the convergence properties of Stieltjes continued fraction approxi- 
mation are needed. 

Proof of Theorem [3l First of all, we observe that (I17p follows by expressing 
Cf(t) via (fl~4"|) in the definition (|15p of F(s) and exchanging the order of 
integration. This can be done for any s for which Res > 0, on account of 
the Tonelli-Fubini theorem for the Stieltjes integral (see [TO]), and it is well 
known that the so defined F(s) is analytic in the same region. 

We come now to the proof of statement [2j We notice that the form of 
series (fl~6P recalls the Stieltjes work [9j on the link between the asymptotic 
series, their expansion in continued fractions and the Stieltjes transform. In 
dealing with this subject he met with the moment problem that now bears his 
name. We recall that, given a sequence of positive numbers c n , the Stieltjes 
moment problem consists in checking whether there exists a measure $ on 
R+ such that 

u n d<S>(u) = Cn . (18) 
The solution exists if and only if 

det(A n ) > and det(A n ) > Vn , (19) 



o 
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where the matrices A n , A n are defined by 



dcf 



Cf) 
Cl 



Cl 
Cn+1 



C/i 
Cn+1 

C2n 



A, 



Cl 
C2 



C2 
C3 



C n +1 C n+ 2 



Cn+1 
Cn+2 

C2n+1 



, (20) 



and if the equality in JT9J) is obtained for A^, then it is obtained for A n , 
with n > k and for A n , with n > k — 1. Furthermore, the solution is unique 
if there exists D > such that 



< D n (2n)! 



(21) 



holds (for a proof with more recent methods, see |12j). We point out that, in 
our case, since Cf(t) is analytic about the origin, there should exist D > 
such that the previous condition is satisfied, in view of the Cauchy-Hadamard 
criterion. 

Coming back to our problem, we collect the results of paper |9j of interest 
to us in the following statement: if the function F(s) admits the represen- 
tation M1\) . then the continued fraction which approximates the asymptotic 
expansion \16\) converges to F(s) for Re s > 0; moreover, the Stieltjes mo- 
ment problem for the sequence c„, is soluble. 

The previous statement is proved this way. In equation ([17|) we put 

z = f s 2 , u = f u) 2 , 3>(u) = f a(yju). Then we have the formal chain of 
equalities 



+oo 



d$(u) 
z + u 



F(s) 



CO 



^ + ^ + ... 

z 2 Z ?J 



(22) 



where again $ is a Borel positive measure on R + . Stieltjes showed that 
the continued fraction which approximates the term in brackets at the r.h.s. 
converges to the integral at the l.h.s, and that the Stieltjes moment problem 
for the sequence c n has a solution for the measure Then it is straightfor- 
ward to check that a/ir solves the associated symmetric Hamburger moment 
problem, which means that 



1 



7T 



w 2n da(u;) 



+oo 



2n+l 



. 



Q.E.D. 



Remark 4. Notice that we have incidentally shown that the coefficients c n 
are such that conditions f)19|) are satisfied at any order. This restriction on 
the values of the coefficients, together with the analogous ones imposed by 
the application of the previous theorem to the dynamical variables f^ n \ for 
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all n, already contains the restriction given by Proposition [TJ However, we 
point out that Proposition Q] holds up to a given order even if the existence 
of beyond such an order is not guaranteed. 

In view of Corollary [H the question of the asymptotic behaviour of the 
time autocorrelations can be restated as a regularity problem for the measure 
<£, which solves Stieltjes moment problem for the c n , or, equivalently, for a 
which solves the associated Hamburger moment problem. An answer to this 
question can be given only if the whole sequence of moments c n is known. 
However, it is also of interest to understand which is the qualitative difference 
between the approximation (at a given order) of an analytic measure and 
that of a measure with at least a singular point. We try here to give an 
answer to the latter question, deferring to Appendix [A] the discussion of 
which is the link between the behaviour of the moments and the regularity 
of the solution (see Proposition [21 [3l 2] therein) . 

We have already pointed out that a rational approximation F2 n (s) of the 
Laplace transform F(s) can be constructed from (|22[) by standard methods 
(see Chapter 2 of Stieltjes memoir [9]), which give 

„ / x P2n{s) <A Pk 

Can (a) t^o s +UJ k 
Such a function, in turn, can be rewritten as the following Stieltjes integral 

F2n{s) = s J 7^7' 

in which the measure a n is piecewise constant and has a jump p^ at the 
points u = uj\. Therefore, in order that the limiting measure a = lima n 
be continuous, it is necessary that, as the approximation order grows, the 
poles become dense and their residues tend to zero. As a consequence, we 
will take as a qualitative indication of a sub-exponential decay the property 
that the rational approximations have at least a pole which remains isolated, 
with a stable residue, as the order grows. 

4 Numerical study of the FPU chain 

Checking whether the spectral measure is analytic, or even whether the cri- 
terion of Proposition [2] of Appendix [A] is satisfied for a concrete problem 
is a major task, because it requires information on an infinite number of 
coefficients. Thus, as a first attempt to implement the ideas of Section [3l 
we computed numerically some coefficients c n for a model which is widely 
studied in the literature and is of great relevance for the foundations of sta- 
tistical mechanics, namely the Fermi-Pasta-Ulam model (FPU, see |14j). It 
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describes a one-dimensional chain of N + 1 particles interacting through non- 
linear springs. The Hamiltonian of such a system, after a suitable rescaling, 
can be written as 




(gj ~ Qj+i) 2 
2 



+ 



a (Qj ~ Qj+if . P (Qj ~ gj+i) 4 



3 4 



where p = (po, ■ ■ ■ ,Pn) an d Q = (Qo, ■ ■ ■ , Qn) are canonically conjugated 
variables, while a and j5 are parameters that control the size of the non- 
linearity. We choose as invariant probability measure the Gibbs one, i.e., 

d/x(p, q) = Z _1 (T) exp(—H(p, q)/T)dpdq, where Z(T) is the partition func- 
tion and T > the temperature, and for the purpose of speeding up the nu- 
merical evaluation of the integrals involving the measure, we consider here a 
chain with one end point fixed and the other one free, i.e, we impose only the 
boundary condition qo = 0, po = (see footnote H] below). We note however 
that a few computations made on the model in which both end-points are 
fixed have shown no significant difference. Concerning the parameters, we 
consider the often studied case in which a = f3 = 1/4. The relevant fact is 
that as T approaches the contribution due to the nonlinear terms becomes 
statistically negligible. 

It is well known that if the nonlinear terms are neglected then the Hamil- 
tonian is integrable, admitting iV normal modes of oscillation, and that the 
motion is quasi periodic so that, for every k = 0, . . . , N — 1, the energy 
of the fc-th mode is a constant of motion (and thus, as obvious, the time 
autocorrelation of such a dynamical variable is constant). Our aim is to 
investigate what happens when the nonlinearity is introduced as a perturba- 
tion, namely, for T positive but close to 0. We are particularly interested in 
the exchange of energy among the normal modes. To this end, a good choice 
could be to consider the fraction of energy localized on the low frequency 
modes, defined, for example, as 



since its time variation is a convenient indicator of the flow of energy from 
low to high frequency modes, and vice versa. As a numerical tool we want to 
compute numerically the coefficients of the series (|10p for this function, since, 
in view of Theorem [3l they give indications on the asymptotic behaviour of 
the autocorrelation. 

As a matter of fact, we notice that the function £ may not represent the 
best choice. For, such a function is strongly correlated with the Hamiltonian, 
so that, in virtue of Theorem [TJ its autocorrelation remains close to that of 
the Hamiltonian, thus remaining far from zero. So, we rather took a suitable 




fc=o 
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modification of £ , namely, the projection of £ , via the Gram-Schmidt or- 
thogonalization process, on the space of the dynamical variables uncorrelated 

with H. Thus we consider the dynamical variable £ = f £ — cov(£, H)H/af^. 

Our calculation proceeds as follows. We extract a sample in phase space 
according to the Gibbs measurd^ and we estimate the L 2 -norm \\£^ || 2 of the 
functions £^ = [£^ n ~ l \ H] generated by £^ = £ taking the mean value on 
our sample. For this purpose, we observe that £^ is a certain function of p 
and q, which can be expressed via Faa di Bruno's formula for the derivatives 
of any order of a composed function. This formula requires the computation 
of some combinatorial coefficients, which needs a smart implementation in 
order to reduce the computation time. We report in Appendix ICl the scheme 
we have followed. 

In order to study the Laplace transform of the time autocorrelation of 
£ , following the procedure of Stieltjes, we approximate it up to order n 
with rational functions, and look for the poles of such an approximating 
function in the complex plane (see chapter 2 of the memoir |!5]). Our aim 
is to compare the behaviour of the time autocorrelation of £ with some 
functions which behave in a somehow known way, by comparing the poles of 
the corresponding rational approximations of the Laplace transforms. The 
first comparison function is the time autocorrelation of a variable which 
should lose correlation in a very short time, i.e., (the orthogonal projection 
of) the kinetic energy K of the first half of the chain's particles. The second 
one will be a given function of time which has an exponential decay (see 
below for its actual form). 

Let us describe the results in some detail. We know that the rational 
approximation of the Laplace transform at order 2n will have 2n complex 
conjugated poles, Sk = ioJk, Sfc = —i^k sav - The inverses of the cj's so found 
for the time autocorrelation of £ are plotted in fig. [1] versus temperature T. 
The pole that shows up at order 1 produces the upper curve in the figure. 
It corresponds to oscillation of a quite large period, growing as 1/T 1 / 2 as T 
decreases to zero. The actual value changes a little when the approximation 
increases, keeping however the same order of magnitude. The plotted data 
correspond to the approximation order n = 4. We have also noted that the 
values found do not significantly change with the number N of particles, 
for N up to 364. The frequencies represented by the lower curves show up 
at successive orders n = 2, 3, 4, respectively, and correspond to oscillations 
of a short period, which, at variance with the first period, are of the order 
of magnitude of the typical periods of the normal modes. They appear to 
remain almost constant with T. 

A further information is provided by the residues of the poles, which 

4 We extract pj and r,- == Qj — qj-i, for j = 1, . . . , N, which are stochastically indepen- 
dent variables for the Gibbs measure with respect to the present Hamiltonian. Here lies 
the great advantage of studying the chain with one end-point free. 
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Figure 1: Inverse of the imaginary part of the four poles of the Laplace 
transform of the time autocorrelation of £ versus temperature, in logarithmic 
scale. Here the number of particles is N = 40. 

represent the amplitudes of the oscillations. The quantities for the four poles 
of fig. [T] are reported in fig. [21 after a normalization defined by setting the 
sum of the residues equal to one. The remarkable fact is that the amplitude 
corresponding to the longest period is the largest one. Even for T = 1 the 
amplitudes of the shorter periods do not exceed one percent of the total. 

It is now of interest to look at the corresponding graphs (fig. I3ati3bp for 
the time autocorrelation of the kinetic energy of half of the chain, K. As 
expected, the time-scales involved are much shorter, as is seen by inspection 
of fig. [3al In any case, an indication of a bigger chaoticity is given by the 
fact that all the residues are of the same order of magnitude and are closer to 
each other, as can be seen in fig. I3bl This means that in the spectrum of the 
time autocorrelation of K at least three frequencies are excited, at variance 
with the spectrum of £ , for which just a frequency is actually excited. 

Then, we compare the behaviour of Cg(t) with that of g(t) = f cosh _1 (6t), 
where b > 0, in order to look for possible differences with respect to an 
exponential decay. Here the constant b is so chosen that the approximation 
of the Laplace transform of such a function at order n = 1 has a pole at 
the same place as the Laplace transform of Cg(t), for a fixed temperature 
T (as we have seen, b ~ 1/T 1//2 ). In fig. 2] the poles and the corresponding 
residues of the Laplace transforms of cosh -1 (6t) and Cg(t) are plotted at 
the successive orders of approximation n = 3,4. 
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Figure 2: Residues at the four poles, whose imaginary parts are reported in 
fig. [1] versus temperature, in logarithmic scale. Here the number of particles 
is N = 40. 
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In both pole of order 10 3 has the largest residue, but it is 

apparent that, while for g{t) the other poles are close to that one, in the case 
of Cg(t) there seems to remain a jump between the first pole and the others. 
This fact is in agreement with the remark, already made by Stieltjes, that 
the poles of the Laplace transform of g{t) become dense on the imaginary 
axis and corresponds with the fact that the spectral measure is in this case 
analytic. On the other hand, the same fact seems to suggest a qualitatively 
different behaviour of Cg(t), and in particular that a(u) is not analytic in 
this case, so that the decay of the time autocorrelation of £ is not exponential. 

As a final remark, we add a few considerations concerning the actual 
difficulty occurring in the calculations. In view of Theorem [3l all poles lie 
on the imaginary axis and all the determinants (|20p are positive. So, we 
are certain that the outcome of our computation ceases to be reliable at the 
order at which we encounter a negative determinant. Since the coefficients 
are determined via a Monte Carlo approximation of an integral, they are 
subject to a numerical error. The difficulty is that increasing the order 
requires also a corresponding increase of the precision, and so a longer and 
longer calculation. This makes the whole procedure unpractical even at not 
too high orders, and for this reason we performed most calculations up to 
order 8 only, thus finding 4 frequencies. 
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Figure 3: Location and residues of the first three poles of the Laplace trans- 
form of the time autocorrelation of kinetic energy K versus temperature, in 
logarithmic scale. Here the number of particles is N = 40. 



5 Conclusions 

We have shown that there is a convergent algorithm which, in principle, 
could be implemented to get the Laplace transform of the time autocor- 
relation of a dynamical variable, and so to get conclusions about the rate 
of its decay as t — > oo. In our opinion, the most relevant aspect is the 
connection we found between the spectral measure and the solution of the 
Stieltjes moment problem, which implies a relation between the smoothness 
property of the measure and the fulfilment of an infinite set of constraints 
(see, as an example, Proposition [3] in Appendix [A]) . This leads us to the 
question, whether smooth measures are in some sense exceptional, namely, 
to ask whether these constraints can still be satisfied when small changes 
in the moments are performed (i.e., by generic small perturbations of the 
Hamiltonian flow). 

Implementing such an algorithm on a computer raises major difficulties 
with respect to both the computational time and the reliability of the results. 
As a matter of fact, the main problem is connected with the calculation of 
large determinants. However, we tried to apply the method to two func- 
tions of interest fro the FPU problem, namely, the energy of a packet of 
low frequency modes and the kinetic energy of half of a FPU chain. For 
the first function, our results seem to support the conjecture that the time 
autocorrelations exhibit a sub-exponential decay, while no definite conclu- 
sions can be drawn for the second quantity. We remark that this is a known 
characteristic for chains of FPU type, namely that quantities related to the 
particles typically exhibit a behaviour coherent with the predictions of Sta- 
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Figure 4: Residues of the poles of the Laplace transforms of Cg(t) and of 
g(t) versus their imaginary parts, in logarithmic scale. The circles refer to 
Cg(t) and the triangles to g(t), while the dashed lines indicate that the order 
of approximation is 3, and the continuous lines indicate order 4. Here the 
number of particles is N = 40 and T = 10~ 5 
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tistical Mechanics, while this often does not happen with quantities related 
to the modes. On the other hand, we should emphasize that the existence of 
functions with a slow decay of the autocorrelations is in itself an obstacle for 
the ergodic behaviour of a dynamical system. For similar results, see, e.g., 

m. 



A Conditions on the regularity of the measure which 
solves the moment problem 

The problem of determining a necessary and sufficient condition in order that 
the solution of an Hamburger moment problem admit a bounded positive 
derivative has been solved by KHz and Krein (see [17J). Unfortunately, their 
criterion cannot be applied to study the derivative itself. Thus, we report 
the criterion of Akhiezer and Krein for checking the existence of a derivative 
for the Hamburger moment problem for the measure a, then we give an 
alternative procedure for studying the derivatives of -^/^'(^Ifl at least in 
the case in which the time autocorrelation is analytic on the whole real axis. 
We convert the problem into an Hausdorff moment problem (i.e., a moment 
problem on a finite interval), in order to get condition on the regularity of 
the successive derivatives. 

To state the result of |17) of interest to us, let us introduce the definition 
of a non-negative definite sequence of moments c n , by requiring that 



Ci+jXiXj > 



»J=1 



for any n and any sequence (X±, . . . ,X n ) of complex numbers jf| Then 
has 



one 



Proposition 2 In order that a'{oj) exist and be bounded, it is necessary 
and sufficient that there exist L > such that the sequence to(L),ti(L), . . . 
defined b$\ 



t k (L) 



(k + l)\L k + l 



2c; 



-L 



k-l 



(fc-iK- s 



kc' 



k-l 



(k 
(k 





-2L 

2K-3 
l)Cfc_2 



o 
o 

-kL 



is non-negative definite. 



Recall that a'(uj) is analytic if and only if is (see Corollary [TJ) . 

6 Notice that this implies that the determinants of A n defined in (|20[) are all non- 



negative, but the converse is not true. 

7 Recall that we have defined c' 2rl = c„ and c' 2n+1 



0. 
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For the study of the successive derivatives, we observe that in paper |13| 
Hausdorff found a solution method for the moment problem on the interval 
[0, 1], which also allows to give a criterion to check whether, at any order, the 
derivatives of the measure exist and are bounded. Such results can be also 
applied to our problem if the time autocorrelation is analytic on the whole 
real axis, because, through an analytic change of variable, one can map [0, 1] 
on R + and obtain the next Proposition [HJ which provides necessary and 
sufficient conditions in order that the n-th derivative of y/u<&'{u) exist and 
be bounded. For its statement, we need to define the auxiliary moments 
and we need to preliminarily define 

def f +00 d$(u) 







z + u 



then we set 



def def 1 d t I - , \ c , , , 

W) = c Hk = co - T , 7TT il -i z ${z) , for k > , 



(k- 1)1 d k 



z 



n=0 ^ n ' 1=0 ^ ' 



(23) 



These quantities are well defined, since it is known that $(z) is analytic in 
z = 1 and the binomial expansion converges]! Then, we define the quantities 



a: 



A^ m = (p+l)(A^ m 1 -Aj i - m 1 _ 1 ) , forA;>0,A ; <m<p-A ; , 
and state the following 

Proposition 3 Let C/(t) be an analytic function of t, for t & R, and let $ 

be the measure which solves the related Stieltjes moment problem. A bounded 
derivative of order n of \pu&'(u) exists if and only there exist > such 
that, for any p>0,0<k<n and k < m < p — k, one has 

(p + l)|Aj )m | <L k . (25) 

Proof. The proof is performed through the Euler transformation w = f ^xj. 
This enables one to define the positive Borel measure ^f(w) on [0,1], by 

def 

ty(w) = 3>(w/(l — w)). The relation between the moments of ^/ on the 
interval [0, 1] and the moments c n of <3? can be obtained by the following 
chain of equalities, for k > 0, 



See the discussion in the proof below. 
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as can be seen by expanding u k = (u + 1 — u) k via the binomial formula. It 
is easy to see that this implies that the coefficients defined by (|23p are 
the moments of \£. 

In order to study the derivatives of a'{uo) = 2^/u<&'(u), we pass to the the 
moments given by the distribution function \fw^'{w). This function exists 
if Cf(t) is analytic on the whole real axis, because in such a case its spectral 
measure a(ui) tends exponentially fast towards its limit for u) — > oo and, in 
consequence, ^ {w) tends towards the same limit as w — )■ 1 at least as fast 
as e ~ K N 1 ~ w _ Moreover, y/w^f'(w) admits a bounded derivative of order n if 
and only if y/u&(u) does. We need, thus, to find the moments of ^/w*$>'(w), 
which are given by 

£ w k+1 / 2 <S>'(w)dw = £ f; + J /2 ) 2 k+1 ' 2 ~ n (w~iy &(w)dw , 



in view of the binomial expansion of \/l/2 + (u> — 1/2). Since this expansion 
converges uniformly for w £ [0, 1], one can exchange the order of summation 
and integration and prove that these moments are indeed the /2fc defined by 
([231). 

As Hausdorff proved, in order that the moment problem for the sequence 
/2fc have a solution ipo on [0, 1], which is almost everywhere differentiable with 
a limited derivative (in our case, it is ip\{w) = y/w'^'{w)) 1 it is necessary and 
sufficient that (p+ l)|Ap m | be bounded, where Ap m are defined by equation 
(|24p . This corresponds to 

A ° /m = CO {^(^^"^h • 

Then, one considers the Hausdorff moment problem in which the measure 
(not necessarily positive) is ipi, and expresses its moments fll through the 
original moments by integrating by parts. So one obtains that the coef- 
ficients Ap m defined by f|24[) play for ip\ the same role which the m play 
for ipo, because, for 1 < m < p — l|^| 

A^, m = (p + l)(A° m - A° m _!) = ( P + X ) jf 1 u; m (l - W r- m d4(») . 

Thus, one can express the condition that ipi admits a limited derivative by 
asking that (p+ l)|Ap m | have an upper bound. The proof is then concluded 
by iterating this way of proceeding. 

Q.E.D. 

9 We neglected here the terms for m = and m = p, since they give indications only 
on the discontinuities of ipi, which can be disposed of, since ipi can be defined at will on 
a set of measure zero. 
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Finally, we add the following Proposition, which can be proved by the 
easy remark that a power series with positive coefficients, which has a radius 
of convergence p, has a singularity for z = p (this is known as Vivanti's 
theorem). 

Proposition 4 In order that the decay of the time autocorrelation of the 
dynamical variable f be not exponentially fast it is sufficient that there exists 
p > such that limsup n _ >+00 ^/c^ < p . 

Remark 5. The condition of Proposition [4] requires only to have an upper 
bound on the c n , whereas those of Propositions [2E] ask for a more detailed 
knowledge on the coefficients c n . We believe that the first one is seldom 
fulfilled (an example in which it is fulfilled is that of a harmonic oscillator). 
For, the n-th derivative of a function usually grows as n!, so that 

lim sup -y/c^" = +oo . 

n— >+oo 

The other requirements enables us to deal with a larger class of variables, 
but, as just said, it has the drawback of needing a very detailed knowledge 
of all coefficients. 



B Proof of Proposition Q] 

Consider the series ([5]), and recall that the remainder has a definite sign. 
Due to the inequalities ([2]) we readily get the bounds 

2n 2n+l 

a° + Y,i-) k 4t 2k >-a° and a° + ^ (-) k a° k t 2k < a° , Vi € M , 

k=l k=l 

which hold true because the remainder is negative in the first case and pos- 
itive in the second one. Setting y = f i 2 , we get the inequalities 

2n 

P n (y) ^ 2a° + ^(-) fc a°/>0, (27) 

k=l 

2n 

Qn(y) = B-)V>0- (28) 

k=0 

i.e., the polynomials P n {y) and Q n (y) are positive semi-definite for y > 0. 
Look now for the real roots of the equations P n {y) = —£ arid Q n (y) = —£, 
with < e < min{ao,ai}- Since P n (y) and Q n {y) ar e polynomials of degree 
In in which In changes of sign occur, Descartes' rule of signs implies that 
such roots, if they exist, must be positive. Thus P n and Q n are positive 
semi-definite on the whole real axis, i.e., the inequalities ([27]) hold true for 
all y G M. 
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Consider now the dynamical variable /W, calculated by recursive application 
of the Poisson bracket with H. By a straightforward application of Theo- 
rem [2] we get that the coefficients of f l are precisely a l n as defined in (fT3j) . 
Since the argument above applies to any dynamical variable, then we get that 
the corresponding polynomials -P^(y) and Q l n (y) are positive semi-definite, 
too, as claimed. 



C Computation of the coefficients needed for the 
Faa di Bruno formula 

Let us recall the Faa di Bruno's formula (see [18j). Let a function /((/(a;)) 
be given, where / and g possess a sufficient number of derivatives. Then one 
has 

where the sum is over all n-tuples of nonnegative integers (k\, . . . , k n ) satis- 
fying the constraint 

1 • k\ + 2 • &2 + 3 • ks + • • • + n ■ k n = n . 

In our case, we apply twice the formula, the first time to find the suc- 
cessive derivatives of the canonical coordinates q and p = q with respect to 
time, the second one to find the derivatives of S (p, q) and K(p, q). Whereas 
the latter application is trivial, we remark that the former is obtained by 
explicitly expressing the force q = p as a function of q and p, then by writing 

d n d n-2 

As the term at the r.h.s. involves only derivatives of q up to order n — 1, the 
scheme can be applied iteratively to obtain all the derivatives of q (and p) 
with respect to time. 

Since the derivatives of / and g are easily calculated in our case, the 
problem is just to have an effective algorithm that produces all the required 
n-tuples (ki, . . . ,k n ). 

For positive integers n, s we introduce the sets 

£n,« = {k £ Ng : ki + 2k 2 + . . . + nk n = s} 

where No is the set of non negative integers. For n = 1 we have /Ci )S = 
{(s)}, i.e., just one element. For n > 1 the following statement is obviously 
true: every element (fci, . . . ,k n ) G JC ns can be written as (k',k n ) where 
k' = (ki, . . . , k n -±) £ K>n-i,s-nk n ■ This is the key of our algorithm. 
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In KL n)S we consider the inverse lexicographic order, right to left. More 
precisely, the ordering is recursively defined as follows: for n = 1 the ordering 
is trivial, since there is just one element; for j, k £ K. n>s with n > 1 we 
say that j < k if either case applies: (i) j n > k n , or (ii) j n — k n and 
(j'l, • • • , jn-i) < (fa, ■■■ , k n -i)- E.g., setting n = 3 and s = 5we get the 
ordered set 

{(0,1,1), (2,0,1), (1,2,0), (3,1,0), (5,0,0)} . 

Remark that the last vector k £ /C„ iS is (s, 0, . . . , 0). 

We come now to stating the algorithm. We use two basic operations, 
namely: (a) find the first vector k £ K. n>s according to the ordering above; 
(b) for a given k £ K, n ,s find the next vector. 

The first vector is easily found by setting 



s 


j fai-l — 


s - nk n 






n — 1 


-71- 





, . . . , ki = s - nk n — . . . — 2&2 • 



For a given k the next vector is found as follows: starting from k<i find the 
first k m which is not zero, so that fc 2 = . . . = fc m _i = and k rn / 0. Then 
replace k m with k m — l and (ki, . . . , k rn -i) with the first vector of K, m -\,k x +m ■ 
The algorithm stops when k = (s, 0, . . . , 0), namely the last vector, because 
the index m can not be found. 
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